# -*- coding: utf-8 -*-
"""
Created on Fri Nov  9 20:57:38 2012

@author: nico
"""
import numpy as np

#from tptsig.wavfile import *
#from tptsig.nmf import *
from scipy import *
from scipy.signal import *
import pylab as plt
from importfilters import *




def ZeroCrossingTDO(filename):
    """
    Recover instantaneous TDO signal by detecting zero crossing
    
    Parameters
    ==========
        Filename of a Gagescope file,
        TODO:
        Sample rate of the datas,
        Number of column of the datafile
        
    Returns
    =======
        Plots of the recovered signal
        Numpy array containing datas
        
    Examples
    ========  
    >>> filename="/Users/Nicolas/Documents/Projets/TDO/post traitement/tracking freq/s120725d.sig"
    >>> temp=loadfile(inputfile,10E6,4)
        
    """
    y=loadfileGage(filename)
    
    # To reduce signals size
    # y=y[:5002]
    
    
    # Detect zero crossing indexes
    # TOTEST: old zc = np.where(np.diff(np.sign(y))=1)[0]
    zc = np.where(np.diff(np.sign(y))>0)[0]
    
    # raw intervals datas
    #ints=zc[range(1,len(zc),1)]-c[range(0,len(zc)-1,1)]
    
    # calculate the zero crossing point
    by=zc+(abs(y[zc]))/(abs(y[zc]-y[zc+1])+0.0)
    
    
    # Calculate intervalls
    bya=by[:len(by)-1]
    byb=by[1:]
    intsb=byb-bya
    
    #average positives and negatives intervals (useless if an entire period is taken)
    #intsb2=(intsb[range(0,len(intsb)-1,2)]+intsb[range(1,len(intsb),2)])
    figure()
    plot(intsb)
    
    
    # Data filtering
    cutoff = 100.
    fs = 1.
    nyq = fs/2.
    filterorder = 5
    
    b,a = filter_design.butter(filterorder,0.1,btype='low') 
    fints=filtfilt(b,a,intsb)
    plt.plot(fints)
    
    return fints


def FFTZeropadTDO():
    (Pxx, freqs, bins, im)=plt.specgram(temp, NFFT=4096, noverlap=1024)
    plt.plot(Pxx[:10,:1000].T)
    y=Pxx[:,10]
    x=freqs
    
    plt.plot(freqs,log(y))
    
    
    _max, _min = peakdetect(y)
    
    print _max, _min
    
    xm = [p[0] for p in _max]
    ym = [p[1] for p in _max]
    xn = [p[0] for p in _min]
    yn = [p[1] for p in _min]
        
    plt.plot(x,y)
    plt.hold(True)
    plt.plot(xm, ym, 'r+')
    plt.plot(xn, yn, 'g+')
        
    # manual spectrogramm
    Nbslices=10000
    # slices
    SlicesDim=(Nbslices,size(y0)/Nbslices)
    y0.resize(SlicesDim)
    
    zeropad=2**12
    
    s1=size(np.fft.rfft(y0[1,:],n=zeropad))
        
    stfft=sp.zeros((s1,Nbslices),dtype=complex,)
        
    for i in range(Nbslices):
        stfft[:,i]=np.fft.rfft(y0[i,:],n=zeropad)
        
    power=log(abs(stfft))**2    
        
    plt.imshow(power[:300,:5000],origin='lower',aspect='auto')
    
    
filename="/Users/Nicolas/Documents/Projets/TDO/post traitement/tracking freq/s120725d.sig"
ZeroCrossingTDO(filename)